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EIGENFUNCTIONS AND THE DIRICHLET PROBLEM FOR THE 
CLASSICAL KIMURA DIFFUSION OPERATOR 

CHARLES L. EPSTEIN* AND JON WILKENINGt 

Abstract. We study the classical Kimura diffusion operator defined on the n-simplex, 

L Klm = Xi(5ij — Xj)d Xi d X j, 

which has important applications in Population Genetics. Because it is a degenerate elliptic operator acting on a 
singular space, special tools are required to analyze and construct solutions to elliptic and parabolic problems defined 
by this operator. The natural boundary value problems are the “regular” problem and the Dirichlet problem. For the 
regular problem, one can only specify the regularity of the solution at the boundary. For the Dirichlet problem, 
one can specify the boundary values, but the solution is then not smooth at the boundary. In this paper we give a 
computationally effective recursive method to construct the eigenfunctions of the regular operator in any dimension, 
and a recursive method to use them to solve the inhomogeneous equation. As noted, the Dirichlet problem does 
not have a regular solution. We give an explicit construction for the leading singular part along the boundary. The 
necessary correction term can then be found using the eigenfunctions of the regular problem. 

Key words. Kimura Diffusion, Population Genetics, Degenerate Elliptic Equations, Dirichlet Problem, Eigen¬ 
functions 

AMS subject classifications. 35J25, 35J70, 33C50, 65N25 

1. Introduction. An n-simplex, E„, is the subset of K n+1 given, in the affine model , 
by the relations: 

(1.1) x\ + • • • + x n+ i = 1 with 0 < Xi for 1 < i < n + 1. 

In population genetics problems, a point (xi,..., x n +i) G E„ is often thought of as repre¬ 
senting the frequencies of n + 1 alleles or types. More generally, E n can be thought of as the 
space of atomic probability measures with n + 1 atoms. Mathematically, E„ is a manifold 
with corners; its boundary is a stratified space made up of simplices of dimensions between 
0 and n — 1. 

The Kimura diffusion operator, which acts on functions of n + 1 variables, 

(1.2) L Klm = ^ Xiffiij - Xj)d Xi d Xj 

l<i,j<n+l 

appears in the infinite population limit of the ( 71 +l)-allele Wright-Fisher model. It represents 
the limit of the random mating term, and actually appears in the infinite population limits of 
many Markov chain models in population genetics, see [12, 30, 16, 5]. The Kimura diffusion 
operator has many remarkable properties, which we employ in our analysis. The properties 
of this operator reflect the geometry of the simplex in much the same way as the standard 
Laplace operator reflects the Euclidean geometry of R”. 

To include the effects of mutation, selection, migration, etc. the operator is modified by 
the addition of a vector field 

n+1 

(1.3) V = ^2bj{x)d Xj , 

i=i 
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which is tangent to E„, and inward pointing along the boundary of the simplex. In appli¬ 
cations to population genetics, the coefficient functions {(.:/;)} are typically polynomials. 
Linear terms usually suffice to model migration and mutational effects, whereas higher order 
terms are needed to model selection, see [30, 33, 39]. 

There are many statistical quantities of interest in population genetics that can be com¬ 
puted by solving boundary value problems of the form 

(1.4) (L Klm +V)u = / in intE„, with u\ bSn = g. 

For example, the probability of a path of the underlying process exiting through a given 
portion of the boundary, or the expected first arrival time at a portion of the boundary, are 
expressible as the solutions of such boundary value problems. Examples of this type can be 
found in [34], and are further discussed below. 

A method for solving some of these problems, at least in principle, is given in [38], 
though it is not very explicit. In this note we introduce a computationally effective method 
for solving high-dimensional, inhomogeneous regular and Dirichlet problems for the Kimura 
operator itself, i.e. with V = 0. For the regular problem, one specifies minimal regularity 
requirements for the solution at the boundary. For the Dirichlet problem, the boundary values 
are specified, but the solution is then not smooth at the boundary. Our method also clarifies 
the precise regularity of the Dirichlet solution in the closed simplex, at least when / and g 
are sufficiently smooth. In addition we show, in somewhat greater generality, how to find the 
eigenfunctions of these operators, which are represented as products of functions of single 
variables, and how to compute the expansion coefficients. 

The operator, L Klm +V, and variants thereof, appear in many classical papers in popula¬ 
tion genetics, see [12, 13, 27, 26, 28, 45], Recently, there has been a resurgence of interest in 
using the Kimura diffusion equation as a forward model for maximum likelihood estimators 
of selection coefficients, demographic models, mutation rates, effective population sizes, etc. 
The evolution of other observable measures of genetic variability such as the allele frequency 
spectrum, or site frequency spectrum, can also be shown to satisfy a variant of the Kimura 
diffusion equation, see [11, 19]. To use this diffusion process as a forward model, one either 
needs to have efficient means for solving the Kimura diffusion equation, see [19,41, 3, 32], or 
one must simulate the underlying stochastic process, [11, 10]. In most previous work where 
the Kimura diffusion equation is solved, the underlying space is 1-dimensional. Even in one 
dimension, many authors employ numerical methods that rely on finite difference approxi¬ 
mations. These are, however, not reliable for imposing the subtle boundary conditions that 
arise with degenerate operators like L Klm , and can in fact lead to errors; see [24], and the 
supplement to [2], By contrast, our approach provides a stable construction, mathematically 
equivalent to a Gram-Schmidt procedure, of bases of eigenfunctions. These can be used to ac¬ 
curately solve both elliptic and parabolic problems, as well as compute approximations to the 
heat kernel itself, which are of central importance in a variety of applications; see [40,41,42]. 

Our construction for the polynomial eigenfunctions of L Klm +V is applicable provided 
that the operator has “constant weights,” see (2.15) below. The case of positive weights 
has been studied extensively in the literature. Kimura [29] and Karlin and McGregor [25] 
used hypergeometric functions to study the two-dimensional case. Littler and Fackerell [35] 
generalized Karlin and McGregor’s approach to higher dimensions using bi-orthogonal poly¬ 
nomial systems. Griffiths [16] showed that the polynomial eigenfunctions of the diffusion 
operator corresponding to a repeated eigenvalue can be orthogonalized and grouped together 
to form reproducing kernel polynomials that appear in the transition function expansion of 
the diffusion process. One way to represent the orthogonal polynomials in the reproducing 
kernels is via a triangular construction of Proriol [36] and Koornwinder [31] for multivari¬ 
ate Jacobi polynomials. Griffiths and Spano also discuss this construction [17], and provide 
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probabilistic connections to multivariate versions of several families of classical orthogonal 
polynomials, including the Jacobi polynomials that arise here [18]. In the present work, we 
present a direct construction of polynomial eigenfunctions for the V = 0 case, so that it is not 
necessary to take limits of the positive weight case as the mutation rates approach zero [16]. 
One novelty that arises is that when V = 0, as functions on the n-simplex, these eigenfunc¬ 
tions do not belong to a single L 2 -space, but each belongs to an /.—space of some stratum of 
the boundary. Their coefficients in the representation of a function in terms of this spanning 
set are computed as inner products on these lower dimensional strata. 

The constructions presented here can serve as the foundation for a perturbative method 
for solving Kimura-type diffusions with a more complicated vector field, also modeling se¬ 
lection. We will return to these and other elaborations of the theory presented here in a 
subsequent publication. 
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2. Some Facts about the Kimura Diffusion Operator. We begin our analysis by re¬ 
viewing some of the remarkable properties of L Klm . A very important fact about L Klm is 
the result of K. Sato [37], which states that if U is a C 2 -function that vanishes on £„, then 
[L Klm C/][e„= 0 as well. Thus we can start with a function U(x i,.. ., x n +i) defined on 
and extend it to be independent of any one of its arguments, say Xj . This uniquely defines a 
function on the projection of the simplex to the hyperplane {Xj = 0}, and vice versa. For 
definiteness we take j = n + 1. This gives a function 

(2.1) u{x i, ...,x n ) = U(x i, • • • ,x n , 1 - (xi H-b x n )) 

defined on a projective model of the simplex: 


(2.2) = {i £ I" : x\ -)-+ x n < 1 with Xi > 0 for i = 1 ,..., n}. 

To compute L Klm U\s n , we can apply L Klm in n-variables to u : 



i 


Here and in the sequel, when using a projective model, we let x n +\ = 1 • -+x n ). The 

projective model is useful for computations, whereas the affine model shows that this operator 
is entirely symmetric under permutations of the variables (x \,..., x n+ \). In particular, it 
makes clear that, from the perspective of L Klm , all the vertices of £„ are “geometrically” 
identical, something that is not evident in the projective model. 

In this note we further investigate some remarkable properties of the operator L Klm , 
which are hinted at in [38, 39]. We first consider a recursive construction of the polynomial 
eigenfunctions of L Klm . After that we show how to solve the inhomogeneous Dirichlet prob¬ 
lem for L Klm , given in (1.4). It follows easily from Sato’s theorem that this problem does not 
generally have a solution in C 2 (E n ), even if / and g are both in C°°. Our method of solution 
exhibits the precise form of the singularities along the various boundary strata. 

2.1. Vector fields. A generalized Kimura diffusion in £ n , with “standard” second order 
part, is a second order differential operator of the form 


n 



(2.4) 
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The vector field is normally required to be inward pointing, which means that 

(2.5) bj( x)r Xi =o > 0 for j = l,...,n, 

and 

n 

(2-6) ^ ~2'b j (x)\ Xl+ ... +Xn = 1 < 0. 

i= i 

We denote the first order terms by 

n 

(2.7) V = Y}i(*)8*r 

3 =1 


In the affine model 

n+1 

(2.8) L = ^2 Xi(Sij - Xj)d Xi d Xj + V, V = ^ bj(x)d Xj , 

l<z,j'<n+l j=l 

with the additional condition that 

n+1 

(2.9) ^ ' bj(a;)|' a ;i4—bxn+i=i = 0. 

j=i 

Considerable generalizations of this class of operators and spaces are introduced in [7]; 
though in the present work we concentrate on the classical cases of model operators on sim- 
plices and positive orthants, R” . In most of the Probability and Population Genetics literature 
a different normalization is employed, namely L Klm = 2 (X+(+' — x j)d Xi d Xi ). In this 
paper we use the normalization more common in mathematical analysis, which omits the 
factor of 1/2. 

The boundary of E n is a stratified space. If K is a boundary face of E„ of codimension 
n — k, then it is again a simplex and is represented, in the affine model, by a subset of 
the form + ■ ■ ■ + Xi k+1 = 1 with Xi t > 0 for l = 1,..., k + 1, which implies that 
Xj 1 = ■ ■ ■ = Xj n _ k = 0. We setl = {*i,... ,ik+i} and 

(2.10) J = {ji,..., jn- k } = {l,...,n+l}\J. 

We denote this boundary face by Kj- 

Every boundary face is a simplex, and the formula for the operator analogous to L Klm is 
the same for any boundary face. For example for Kx, the sum in (1.2) is simply restricted to 
the variables { x ^,..., Xi k+1 }. We denote this operator by 

(2.11) Lj im = ^2 Xi(5ij - Xj)d Xi d Xj . 

ijei 

Let v be a C 2 -function on Kx and let v denote any C 2 -extension of v to the ambient M" 4 1 . 
K. Sato, in fact, proved, (see [37]) that 

( 2 . 12 ) Lf m v={L Kim v)\ Kl ■ 

Thus, the restriction or extension of L Klm to boundary strata or higher-dimensional simplices 
is canonical, similar to the connection between the projective and affine models discussed at 
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the beginning of this section. The extension property makes it more natural to label the “last” 
variable as x n+ \ rather than xq, since it may not actually be the last. 

In the affine representation, a vector field V (satisfying (2.9)) is tangent to Kj provided 

that 

(2.13) Vxj\ Xj= o = 0, (jeJ). 

Condition (2.9) requires that V be tangent to £„ itself, i.e. V (xi + • • • + x n +i) = 0. For 
vector fields, the analogues of Sato’s results are obvious: if U vanishes on £„ then VU = 0, 
so extending a function U defined on £„ to be constant in any coordinate direction leads to 
an unambiguous value of VU\s n ; and, defining 

(2.14) Vi = Y J h(x)d Xi , 

ie i 

we see that if V is tangent to Ki then Vzv = (Vv) \k x , where v and v are defined as in 
(2.12). 

2.2. Constant weights and the Hilbert space setting. In addition to L Klm there are 
other classes of “special” Kimura diffusion operators. Classically one singles out operators 
with constant “weights.” (This terminology is introduced in [8].) These operators have the 
property that the functions bj[x) (or bj(x)) are linear and 

(2.15) ( Vxj)\ Xj = 0 = bj(x)\ Xj =o= bj, 

a constant. In the projective model, one may readily show that such a vector field takes the 
form 

(2.16) V b = EU {bj-B Xj )d Xj , 

where 

(2.17) b = (bi,, b n+ i) and B = bi H- \-b n+ 1 . 

Note that b n+ 1 enters (2.16) through B. We define 

(2.18) L K im = L Kim + V b . 

These operators are special because they are self-adjoint with respect to the L 2 inner product 
on £ n defined by the following measure, which (up to normalization) has the Dirichlet density 
representing the stationary distribution of the underlying Markov process (see e.g. [16, 17, 
18]) 

n+1 

(2.19) d/ib(x) = w b (x) dx i • • • dx n , w b {x) = ]^[ x- . 

j =i 

Indeed, in this case L^" 11 in (2.8) may be written 

T Kim _ - d Xj )[w b XiXj{d Xi - d Xj )u] 

L b U ~ 2^ 

i<j 


( 2 . 20 ) 


w b 
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where the sum is over all pairs i,j £ { 1 ,..., n + 1 } with i < j , and, in the projective model 

^ ^ n+1 

(2.21) (L «) = (u, Lj+u) = - 

i<j 

Here u, v are independent of x„+i and agree with u, v on In deriving (2.21), it was 
assumed that WbXiXj[(d Xi — d Xj )u]v = 0 and WbXjXj[{d Xi — d Xj )v\u = 0 on the faces 

{xi = 0 } and {x :) = 0 }, so care must be taken in defining the domain of Lj)" n when working 
in the Hilbert space setting with b = 0, see [38]. Note that x„+i is treated as an independent 
variable in these formulas when computing partial derivatives, but x n+ i = 1 — ( x\ + ■ ■ ■+x n ) 
when evaluating integrals over £„. A useful variant of (2.20) is 

\ ' dxi(.^bXiX n -\-iO X j j U') 

r-f w b 

i=i 

where x n +i is now treated as a dependent variable, i.e. <9+i / dx, = — 1. If all the bj are 
positive, then dfib{x) has finite total mass. In this paper we are primarily interested in the 
case where all of the bj vanish, in which case the volume of £„ is infinite. 

2.3. The Dirichlet problem and alternative function spaces. In section 6 we present 
a method for solving the inhomogeneous Dirichlet problem in (1.4). These results have ex¬ 
tensions to the case of Kimura diffusions with constant weights, though Dirichlet boundary 
conditions are only appropriate on faces {xj = 0 } for which bj < 1 . For simplicity, we 
focus our attention here on the case when all the weights are zero, which is already of central 
importance in applications. We also assume that / and g in (1.4) are sufficiently smooth. A 
different analysis of this problem, employing blow-ups, appears in [20, 21, 22]. The blow-up 
approach gives a much less explicit description of the singularities that arise when the data is 
smooth on the simplex itself, but allows for considerably more singular data. 

We use the notation and definitions of various function spaces given in the monograph [7], 
The principal symbol of the Kimura diffusion operator, 

(2.23) P Kim (0 = ^ Xitfy - XjMj, 

i,j 

defines the dual of the intrinsic metric on the simplex. This corresponding metric is singular 
along the boundary and incomplete, i.e., the boundary is at a finite distance from interior 
points. The distance function on E n defined by this metric is equivalent to 

n+1 

(2-24) pi(x, y) = I V*j - VVJ\ ■ 

3 =1 

We also use two scales of anisotropic Holder spaces, C^p(£„) and C^p +7 (£ ra ), k £ No, 7 £ 
(0,1), introduced in [7], with respect to which the operator L Klm has optimal mapping proper¬ 
ties. These spaces are defined with respect to the intrinsic metric. For example / £ C^}p(£ n ) 
if / € C°(£„) and 

rn \f(x)~f(y)\ ^ 

(2.25) L/Jo, 7 ,wf = sup - 7 -r-— < OO. 

/6E„ Pi [X 1 y ) 

If A < 0 and / £ p(£ n ), then the elliptic equation, (L Klm —A )u = /, has a unique solu¬ 

tion, u £ C^)p +7 (£„), indicating that this is indeed the correct notion of “elliptic regularity” 
for operators of this general type. 


( 2 . 22 ) 


■ = £ 
%<j 




W b 


%iXj\(d X i )^] dxj')v\d/Jib(x')- 
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3. A 1-d Example. We begin by considering the Dirichlet problem in the ld-case. 
These results are well known, and serve as motivation for our subsequent development. Sup¬ 
pose that we would like to find the solution to 

(3.1) x(l — x)d x u = f with u(0) = go , m(1) = gi■ 

We can write u = uq + (1 — x)go + xgi, where uq vanishes at the boundary of [0,1]. It is 
apparent that uq cannot be C 2 up to the boundary unless /(0) = /(1) = 0. In fact, 


(3.2) moO) = x logxf(0) + (1 — a;) log(l — x)f(\) + uo{x), 

where 

(3.3) x(l - x)d x u 0 = f(x) - [(1 - x)f( 0) + xf( 1)] = / (1) with u 0 (0) = u 0 (l) = 0; 

the right hand side, is as smooth as /, and vanishes at 0 and 1. As was shown in [ 6 , 7], 
this equation has a unique solution, with optimal smoothness measured in the anisotropic 
Holder spaces, C fc ’ 2 + 7 ([0,1]). In particular, if / e C°°([0,1]) then so is uq. 

The eigenfunctions of x(l — x)0 2 r that vanish at the boundary are polynomials of the 
form x(l — x)pm(x), where p m is the polynomial of degree m > 0 that satisfies the equation 

(3.4) [x(l - x)dl + 2(1 - 2x)d x - 2 ]p m = Xmp m . 

For later reference, the left-hand side may also be written [I.!) 11 " —2 \p m . By inspecting the 
action on Vd/Va- 1 , where Vd is the space of polynomials of degree at most d, we see that 
A d = —(d+l)(d + 2). Since the eigenfunctions are orthogonal with respect to dpo{x), these 
polynomials are orthogonal with respect to the inner product 

(3.5) ( p,q)=l p(x)q{x)x a (l - x^dx, 

Jo 

with a = fi = 1. Thus, they are multiples of the corresponding Jacobi polynomials [43, 14], 
Pd{x) oc P^ 1 ' 1 \2x — 1). If we choose the normalization \\pd\\L 2 {[o,i]-,d^) = 1- the result 
would be the same as performing Gram-Schmidt on the monomials { 1. x, x 1 .... }. This 
yields the 3-term recurrence 

Po{x) = Vl/ 70 , Vhpi(x) = (x- ao)po(x) 

(3.6) - — 

\/b m+1 p m +i{x) = (x - a m )pm(x ) - y/bmPm- 1 (%), (m > 1), 


where 70 = 1/6, a m = 1/2 and b m = 4 [ 4 (^+i)a!_ 1 ] hi this case. Solving (3.3) for uq now 
boils down to representing in the eigenbases: 

OO 

= X! ( C m/X m )x(l - X)p m (x). 

m =0 


fW(x) * 

(3.7) f(x) = — -- = ^ CmPm(x) 

- x ) ^0 


u 0 (a 


The simplest way to obtain an approximation of f(x) in Vn- 1 is to evaluate f^\x) at the 
zeros Xj of j)\[(x), and to compute the coefficients using Gauss-Lobatto quadrature: 
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Here we used the fact that / W(:r) vanishes at xq = 0 and x^+i = 1. The abscissas Xj and 
weights ujj are easily found using a variant of the Golub-Welsch algorithm [14, 15]. Further 
details of our numerical implementation will be given elsewhere. 

We can estimate the size of the coefficients {c m } in terms of the smoothness of /W, 
using the facts that 


(3.9) 


Cm = / f (1) {x)p m (x)dx 

Jo 

[ L 2 ,m p m f {1) dx = [ p m [L Kim +2]f^dx. 
Jo Jo 


and 


The second formula is a special case of (4.1) below. If /W is in C 2l (\ 0,1]), then we can 
iterate the integration by parts formula to obtain: 


(310) Cm “ [m(m. + 3)]* l Pm ^ LKlm + 2 } l f (1)dx - 

As we show below, there is a constant C, independent of m so that 


(3.11) Ibmlk- < C\J m(m + 3)||p m || L 2 ([0il];dA12) = Cy/m(m + 3)|, 


and therefore, there is a constant C so that 

.|| [L Kim +2] t /( l)|| io 


(3.12) 


\c m \ < C- 


U-i 


[m(m + 3)] 

In this instance it is easy to see that 

(3.13) ||[L Kiiri +2] i / (1) || i oo < C;||[/|| C 2 


This is almost a spectral estimate, but we lose one order of decay, in part because we are 
estimating in the L°°-norm, and in part because there is an implicit division of / 1 ' 1 ’ by x(l—x) 
in the formula for c m . 

We can also use the L 2 -norm to estimate these coefficients via 


(3.14) 


< C‘ 




[m(m + 3)] ; 


While the denominator is now [ 777 ,( 771 + 3 )]* the numerator implicitly involves the L 2 -derivative 
of |[L Klm +2]*/( 1 ^ | 2 at the boundary of [0,1], 

Remarkably, very similar approaches work to find the eigenfunctions of L Klm and solve 
the Dirichlet problem in any dimension. This is explained in the following two sections. In 
Section 4 we give a novel construction for the eigenfunctions of the neutral Kimura diffusion 
on the 77 ,-simplex, which highlights their vanishing properties on subsets of the boundary. In 
Section 6 we show how to solve the Dirichlet problem on an n-simplex, with arbitrary smooth 
data. 

4. The Polynomial Eigenfunctions of L Klm . In this section we give a hierarchical 
method of constructing the eigenfunctions of L Klm , with considerable control over their van¬ 
ishing properties on bT, n . As before, we let Vd denote polynomials of degree at most d: the 
variables involved will be clear from the context. 
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Our results are based upon a formula, which follows easily from a similar calculation in 
the work of Shimakura, see Section 7 of [38], As noted above, we let L^ lm = L Klm +14, 
where 14 has linear coefficients and assigns constant weights to the hypersurface components 
of We then let I = (*i < *2 < • • • < ik+ 1 } C {1,..., n + 1}. Shimakura’s work 
implies the following formula 


(4.1) 

L f' m {w x ip) = w T [L£ im -ki] ip 

Where 

k+1 

(4.2) 

Wi{x) = H ; 


7=i 


(4.3) 


and 


b 'i 


2 — bj if j G {ii,...,i fc+ i}, 
bj if j i {*i, - - -+fc+i}; 


(4.4) k x = 5^(1 - bj ) fc + b 3 ■ 

/ \ m j 

If k = n and b = 0, then b' = 2 = (2,..., 2). Equation (4.1) has a wide range of applica¬ 
tions, and is especially useful in cases where some of the {4} are zero. This formula gives a 
very potent method to construct eigenfunctions that have a simple form and vanish on certain 
parts of the boundary. 

Recall that the C°-graph closure of L Klm acting on C 3 (E n ) is what is called the “regular 
operator” in [7], which is the “backward” Kolmogorov operator in applications to Population 
Genetics. The eigenfunctions that we construct are polynomials and hence in the domain of 
the regular operator. Indeed, we will show that for each natural number d , the eigenfunc¬ 
tions of degree less than or equal to d actually span Vd • Hence this is the complete set of 
eigenfunctions for the regular operator. 

As functions on the n-simplex, these eigenfunctions do not belong to a single L 2 -space, 
but each belongs to an i 2 -space of some stratum of the boundary. Their coefficients in the 
representation of a function in terms of this spanning set are computed as inner products on 
these lower dimensional strata. 


4.1. Hierarchy of Polynomial Eigenfunctions. The simplest polynomial eigenfunc¬ 
tions of L Klm are the functions {xi ,..., x n + 1 }, which are null vectors. Each of these eigen¬ 
functions vanishes on a codimension 1 boundary face of E n . Of course, a constant func¬ 
tion is also a null-vector for L Klm , but it already belongs to the span of the others since 

X\-\ -b x n +i = 1. To fit with the pattern below, one can write these functions as x^ip{x), 

where 'il'(x) = 1 spans the space of constant functions determined by their value at vertex i-\ 
(where x ^ = 1). 

Next, for each distinct pair 1 < i \ < *2 < n + 1, we look for eigenfunctions of the form 
Xi 1 Xi 2 ip{x). If i/>(x) is only a function of (or Xi 2 ) and satisfies [L 2 lm —2 ]ip = Xif), i.e. 

(4.5) 3+ (1 - )dl H ip + 2(1 - 2a+ )d Xii %p-2ip = A ip, 

then, by (4.1), Xi 1 Xi 2 ip(xi 1 ) is also an eigenfunction of the original operator L Klm with the 
same eigenvalue. Note that if {j 1 ,..., j n -i} = {1,..., n + 1} \ {i\, * 2 }, then we are solving 
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on the edge 

(4.6) x :n =■■■ = x Jn _ 1 = 0. 

Hence the eigenfunctions Xi 1 Xi 2 ip(xi 1 ) vanish on the boundary of this edge. Also note that if 
a function ip depends only on a subset of the coordinates, then L^ lm ip also depends only on 
the same subset of the coordinates. In particular, (3.4) and (4.5) agree. 

With these observations, working in the different projective models, we can construct all 
the polynomial eigenfunctions of L Klm . These take the form 

x il * x ik+i^P{ x ii i ■ • * ; •t'ifc ); 

for various choices of indices {ii, ..., ik+i}- Here ip is an eigenfunction of the operator 

(4.7) Lf™ = l*T\ Kx 

acting on a /.'-simplex, and the variables (x,, ,..., Xi k+1 ) range over the face defined by the 
equations 

(4.8) x n =■■■ = x jn _ k = 0, 

using the notation of (2.10). The weights for this operator are all equal to 2, and therefore we 
denote it by . A simple calculation shows that, if \I\ = k + 1, then 

(4.9) L^[x^---xT k k ) =-(\m\ 2 + (2k+l)\m\)(x™'---x™*) mod 

where \ fh\ = mi + • • • + nik- (A more general formula is given in (5.9) below.) Thus, Lj™ 
leaves the subspaces Vd invariant; the dth eigenvalue is [— d 2 — (2 k + 1 )d\ for d > 0; and its 
multiplicity is equal to ( t_Y 1 ), the dimension of Vd/Vd-i- Moreover, applying the Gram- 
Schmidt procedure to the set of monomials in the variables x tl , x lk , ordered by degree 
(but arbitrarily ordered within the set of monomials of the same degree), will lead to an or¬ 
thonormal set of eigenfunctions of L^™. As a result, the eigenfunctions ipx,m ( x h , • • •, %i k ) 
of this operator are multivariate orthogonal polynomials with respect to the measure 

(4.10) d/kr, 2 = (n-i 1 x ij) dx h ‘ ‘ ‘ dx i k ■ 

The corresponding eigenfunctions Wxipi,m of Lj"" are orthogonal with respect to dfi o when 
k = n and I = {1,..., n + 1}, but are not normalizable (i.e. do not belong to L 2 (E n ; d/io)) 
when k < n. Nevertheless, they are still eigenfunctions algebraically, with eigenvalue 


(4.11) Ax i? r = -|m| 2 - (2k + l)|m| - k{k + 1), (fc = \l\ - 1) 

and play an essential role in solving L Klm u = f below, where we use them to adjust / to 
zero on the boundaries, just as was done in (3.2) and (3.3) in 1-d. 

From this observation, it is not difficult to demonstrate the completeness, as a Schauder 
basis for C°(E n ), of the eigenfunctions obtained in this way. Let I = {zi,... ,ik+i\ C 
{1,..., n + 1} be a set of indices, and K'x C bE ra , the corresponding boundary face. We let 
Vi C Cf-x 1 ,;,,..., x lk \ denote the ideal of polynomials defined on Kj that vanish on hKx- It 
is easy to see that 

Vi x i k • • • ,..., x i k ] • 


(4.12) 
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Formula (4.1) shows that this ideal is invariant under L^ 1 "' . The operator is self adjoint 
with respect to the measure ■ ■ ■ Xi k+1 dxi 1 ■ ■ ■ dxi k , and it too preserves Vd , for any d. 
Thus the polynomial eigenfunctions of Lj™ span ,..., Xi k \, and therefore x^ ■ ■ ■ Xi k+1 
times these eigenfunctions spans Vi- 

We next observe that an eigenfunction of this type must vanish on every other fc-simplex 
in the boundary of E„. This is because one of the functions {a;^. : j = 1,..., fc + 1} must 
appear as a defining function for any other fc-simplex. Using this observation, we see that 
any polynomial / that vanishes on every fc-simplex except I\x can be expanded in these 
eigenfunctions. This suggests a simple recursive method for finding the regular solution of 
an equation of the form 

(4.13) L Kim u = f. 

As discussed below, the regular solution is required to belong to an anisotropic Holder space 
in the closed simplex rather than to take on specified boundary values. Even if / is a polyno¬ 
mial and g = 0, the solution to the Dirichlet problem (1.4) is generally not differentiable up 
to the boundary of S„. 

4.2. The Regular Solution of L Kim u = f. As suggested by our construction of the 
eigenfunctions, the regular solution to L Klm u = / is found recursively by working one 
boundary stratum at a time. The fc-skeleton of E„ is defined to be the union of fc-simplices 
in bYi n . We denote this subset by E(j. It is connected for fc > 0, but Ejj \ E^ _1 is a disjoint 
union of open fc-simplices with boundaries contained in Ejj -1 . The regular solution takes the 
form 

(4.14) u = u\ H-+ u n , 

where the term Uk is found by using the eigenfunctions constructed above to solve an auxiliary 
inhomogeneous Dirichlet problem on Ejj \ Ejj -1 . 

If u £ C 2 , then L Klm u vanishes at each vertex. Hence, we start by assuming that / 
is a polynomial that vanishes on each of the vertices of E„. This assumption is dropped in 
Section 6, where imposing inhomogeneous Dirichlet boundary conditions on 6E„ inevitably 
introduces singularities anyway. We define ii \ as the solution of the equation L Klm Mif S l = 
/ [si. As / vanishes on the vertices, which are the boundaries of the components of Ej ( \ 
E°, we can use the eigenfunctions Xi 1 Xi 2 ip(xi 1 ) constructed above to solve this equation on 
each component of Ejj \ E° independently of the others. Note that, using the eigenfunction 
representation, the function u± extends canonically to the entire n-simplex. 

We next solve 

(4.15) L Kim M2 r S 2= f\^ n -L Kim Ul t^ . 

The right hand side vanishes on Ej t , so we can use the eigenfunctions wx'£x,m with X = 
{*i, * 2 j * 3 } to independently solve this equation on each connected component of Ejj \ Ej ( . 
Recursively, we assume that we have found u \,..., Uk- 1 , so that 

(4.16) f-L Kim {u 1 + ---+u k . 1 ) 
vanishes on the (fc — l)-skeleton, and then solve 

(4.17) L Klm Ufc|" S fc = /— L Klm (zti + • • • + . 

Using the eigenfunctions of the form a • • • Xi k+1 ip(xi 1 ,..., Xi k ), we can solve the relevant 
Dirichlet problems independently on each component of E(j \ E(j -1 . 
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The process terminates when we reach the interior of the n-simplex, where we solve the 
problem 

(4.18) L Kim u n h n =fh n -L Kim ( Ul + . 

Once again the right hand side vanishes on the entire boundary of E„ and we can solve this 
equation using eigenfunctions of the form x\ - ■ ■ x n +iif>(x). Since this can be done for any 
polynomial that vanishes on the vertices, it demonstrates that the eigenfunctions constructed 
above, including the nullspace, are in fact a complete set of eigenfunctions for the operator 
L Klm , acting on polynomial functions defined on E n . Since these functions are dense in 
C 2 (£„) it follows easily that this is also a complete set of eigenfunctions for the graph closure 
of L Klm acting on C° (E n ). Altogether we have proved the following result: 

THEOREM 4.1. The regular operator L Klm acting on functions defined on E n has a 
complete set of eigenfunctions of the form 

n 

(4.19) £(L Klm ) = {xi, ..., x n +i} U (J (J x h •••£i fc+1 £( L z,2 1 )- 

fc=l I={l<ii<---<ifc+i<n+l} 

Here £ (L^™) denotes the eigenfunctions of the operator Lj ™, which are polynomials in the 
variables { x j x ,..., Xi k }. 

REMARK 4.2. Another useful consequence of the hierarchical description of the eigen¬ 
functions of L Klm is that it dramatically reduces the work needed to find the eigenfunctions 
o/L Klm acting on E n+ i, once they are known for £„. Indeed, up to choosing subsets of co¬ 
ordinates (xi,..., x n +i), all that is required is to find the “new” eigenfunctions, which are 
of the form 

xi ■ ■ ■ X n+ 2 t/j(Xi, ... ,x n+ i). 


REMARK 4.3. At stage k of the recursive algorithm above, we are given a function f that 
vanishes on the (k — T)-skeleton o/E„ and wish to find Uk such that L Klm Mfc |' s »d = f\s k - 
For simplicity of notation, we have absorbed L Klm (rti + • • • Uk-i) into f in (4.16). Since 
the faces of the k-skeleton decouple, we may assume f = 0 on all faces Kx except one of 
them, which we take to be I = {1,..., k + 1}. Recall from the discussion after (4.9) that the 
relevant eigenfunctions o/L^ lm are of the form (x\ ■ ■ ■ Xk+i)'fix,rh{x), with k 

a set of multivariate orthogonal polynomials with respect to dpx.i hi (4.10) on S; c . Thus, the 
generalization of (3.7) to k dimensions is 
(4.20) 


/(*> = 

«£l * ’ ’ *£/c+1 


mGNg 


Uk{x) 


E 


E-Xi • • • Xfc+1 fix,ih[x). 
A-fh 


In [7] it is shown that for k a non-negative integer, and 7 £ (0,1), the operator L Klm 
maps the anisotropic Holder space C^rp +7 (£„) to C^p(£„). If we let C^p denote the closed 
subspaces of functions vanishing at the vertices of E„, then the inverse 

(4.21) (L Kim ) -1 : (£„) —> C^ 2 F +7 (E n ) 

is shown to be a bounded operator. Let CV 7 denote the bound on this operator. If we approx¬ 
imate / € C^p(£ n ) by a polynomial, p, then 

IKL 17 ™)^ 1 / - (L Klm ) _1 p||wF,fc, 2 + 7 < Co,7l|/ - P||wF,fc, 7 - 


(4.22) 
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The classical (/<:. 7 )-H 6 lder norms dominate the WF-Holder norms, and C C C k . It is 
well known that the accuracy of the best degree cl polynomial approximant only depends on 
the classical Holder smoothness of the data, thus demonstrating the efficacy of our method 
for solving the equation L Klm u = f for any reasonably smooth data / that vanishes at the 
vertices of 

One can generalize the construction of polynomial eigenfunctions to any Kimura oper¬ 
ator, L^ lm , with constant weights. If all of the weights are positive, the measure d/j.^ has 
finite mass and the polynomial eigenfunctions of the regular operator can be chosen to be or¬ 
thonormal with respect to this measure. This case is, in many ways, the easiest to deal with as 
there is a single ambient Hilbert space containing all the eigenfunctions. It has been studied 
previously by solving hypergeometric equations [29, 25], by computing bi-orthogonal sys¬ 
tems of polynomials [25, 35], or by computing orthogonal polynomials on weighted Hilbert 
spaces and grouping them into reproducing kernels [16, 17, 18]. The present work extends 
this third approach to the case of 0 weights, and, for reasons explained below, does not take 
the final step of forming the reproducing kernels. In the construction of Section 4.1 above, we 
look for eigenfunctions of the form X; Ll ■ ■ ■ Xi k+1 ip(xi 1 ,..., Xi k ) with cp an eigenfunction of 
the auxiliary operator L^, lm . If all the weights are positive, the leading factor of x it ■ ■ ■ x lk ( , 
is dropped, k is set to n, and the eigenfunctions of kj)"" are computed directly, without re¬ 
course to an auxiliary operator. While eigenfunctions with respect to a subset of the variables 
are again eigenfunctions of L^ lm , they are already present in the construction at level n, and 
do not have to be dealt with recursively along the stratification of <9E„. Computing these 
eigenfunctions works the same in both cases, and is mathematically equivalent to applying 
the Gram-Schmidt method to the monomials ordered by total degree. This is explained in the 
following section. 

When some weights are zero and some positive, Shimakura’s formula leads to a partial 
hierarchical structure for the eigenfunctions, similar to that described here. This is discussed, 
to some extent, in [38, 39]. As shown in [38] it is possible to specify homogeneous Dirichlet 
data on any boundary face where the weight h < 1 . In general the solution is not smooth 
along this face, but has a leading singularity of the form x\ h '. We will return to these ques¬ 
tions in a subsequent publication. 

5. Numerical Construction of Eigenfunctions. We now give a practical method to 
construct the eigenfunctions described above. By a change of variables these eigenfunctions 
can be represented as products of functions of single variables, which are themselves eigen¬ 
functions of differential operators. Since it is no more difficult, we consider the general case 
of a Kimura diffusion on £„ with constant weights. 

In 2-d, it was observed by Proriol [36] that orthogonal polynomials on a triangle can be 
represented as tensor products of 1-d Jacobi polynomials. See Koornwinder [31] for further 
background and more complicated 2-d geometries. A rather different approach using hy¬ 
pergeometric functions was developed in 2-d by Kimura [29] and Karlin and McGregor [25]. 
This approach was extended to the simplex or n- sphere in the context of the Laplace-Beltrami 
equation by Kalnins, Miller and Tratnik [23]. Littler and Fackerell [35] found solutions of 
Kimura diffusion using expansions in bi-orthogonal polynomials. Griffiths [16] found rep¬ 
resentations for the transition function for this diffusion process using reproducing kernels 
expressed in terms of multivariate Jacobi polynomials; see also [17, 18]. The method we 
derive below is equivalent to Griffiths’ approach, though we avoid computing reproducing 
kernels as it is more efficient and numerically stable (at the expense of symmetry) to leave 
the eigenfunctions separated. Wingate and Taylor [44] showed how to generalize the Proriol 
construction to the A'-dimensional simplex in the case of a uniform weight. We adapt their 
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method to allow more general Jacobi weights of the form 


(5.1) {p,q)= [ p{x)q{x)x^ 1 dxi---dx k . 

The case of interest here is ay = &' — 1 = 1, but treating the general case is no more 
complicated, and can be used to study Kimura diffusion with non-zero weights b , as shown 
below. 

The idea is to map the unit cube to the simplex in such a way that the desired orthogonal 
polynomials on the simplex pull back to have tensor product form on the cube. One way to 
do this, which differs somewhat from the choice made by Wingate and Taylor, is with the 
change of variables x = T[X) given implicitly by 


(5.2) xi = Xi, X 2 = {l — xi)X 2 , ■■■ Xk = {1 - xi- Xk-i)X k , 

which solves to Xj = (nti(i - Xi))Xj and x k +i =l-xi - x k = n-=i(i -Xi). 

These are blow-ups similar to the changes of variables that are used in [20, 21, 22], 

The Jacobian matrix DT(X) is lower triangular, so its determinant is easy to compute: 

(5.3) J= | det DT\ = n- =1 (l - Xi) k - j - 
The inner product (5.1) may now be written 

P k 

(5.4) (p, q)= (p O T)(q O T) ( [] [.Xf (1 - Xj ) {k ~ j)+T -^ “<]) dX 1 ■ • • dX k . 

J[o,i] fc \ =1 ' 

Next we note that if p is a polynomial (in one variable) of degree d, then 

j —i / \ 

(5.5) p{X :j ) [](1 - X i ) d =p( 3 1(1 - si- Xj-xY 

V 1 X 1 X j — 1/ 


is a polynomial of degree d in the variables xj, • • • , Xj. We can therefore define 'Wrri(x) via 


(5.6) 


where 

(5.7) 


3~ 1 


Vte o T(X) = n Qfij’^HXj) JJ(1- Xi) m * J 
j =1 V »=1 / 


3 =1 


k 

Oijm = Ctfc+1 + ( a i + + 1). 

i=j +1 


Here {Q™ ,/3) }“ =0 are orthogonal polynomials on [0,1] in the inner product (3.5), normalized 
to have unit length. Substitution into (5.4) gives 


'tprn') — 



(Xj)xy (i - 


(ai+Trii+m'+l) 
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which, proceeding from j = k to j = 1 , may be seen to equal J|^=i ■ We note that 

iprh{x) is a polynomial of degree d = \fh\ = rn .i■ We order the multi-indices first by 
degree {in < m' if d < d') and then lexicographically from right to left (if d = d! then 
to < to/ if TOfc < m' k , or if m k = rn' k and m k -1 < Tn' k -i' etc., so that (d, 0 , ..., 0 ) is 
smallest). In this ordering, the ipm{x) are the same as one would get from applying Gram- 
Schmidt to the monomials x m in the same order. Hence, all polynomials in k variables are 
accounted for by this construction. For example, when k = 3 and oi\ = «2 = £*3 = 0:4 = 1, 
the first 10 orthogonal polynomials are 

^000 = 12^35, V’ioo = 12^105(4a;- 1 ), ^010 = 12 ^ 210 ( 3 ; + 3y - 1 ), 

ipooi = 36^70(-l + x + y + 2 z), V ’200 = 24^55(1 - 9* + 15a; 2 ), 

V’no = 6V2310(-1 + 5a ;)(—1 + x + 3 y), ip 0 2 0 = 6\/330(3 + 3a ; 2 - 21 y + 28 y 2 + 3*(-2 + 7 y)), 
ipioi = 18\/770(—1 + 5x ){—1 + x + y + 2 z), ^011 = 30%/462(—1 + * + 4y)(-l + x + y + 2 z), 
V ’002 = 24\/ll55(l + x 2 + y 2 — 5z + 5 z 2 + y {—2 + 5 z) + x {—2 + 2 y + 5 z)), 

which span the same space as 1, x, y, z, x 2 , xy , y 2 , xz, yz, z 2 and are pairwise orthogonal. In 
numerical computations, it is best to leave them in the form (5.6) since representation in the 
monomial basis is both expensive and numerically unstable when the degree is large. 

The are proportional to Jacobi polynomials, Qm’^\x) oc Pm' a \2x — 1), but 

are more easily computed directly from a 3-term recurrence. Explicitly, p m {x) = Qm’ P \x) 
satisfies the recurrence (3.6) with 70 = and 

_ 1 /_ (a + P){a - /?) _\ m= 0 a + 1 

am 2 ^ {a + /3 + 2m + 2 ) (a + /3 + 2m) J a + f3- f2’ 

_ m{a + m){/3 + m)(a + /3 + to) TO = 1 (a + 1) (/? + 1 ) 

(a + ft + 2to ) 2 [(a + f3 + 2m ) 2 - 1] (a + [3 + 2 ) 2 (a + /3 + 3) ’ 

These formulas can be derived from the well-known recurrence relations of the Jacobi polyno¬ 
mials [14, 1], modified so that | p rn \j = 1. Finally, now that we have specified the orthonormal 
basis ip^,(x), the coefficients 

(5.8) c A = f(x)i/jYn(x)xi 1 ■ • -x ^ 1 dxt ■ ■ -dx k = / f{x)i/jtf l (x)dx 1 ---dxk 

in (4.20) can be computed by Gauss-Jacobi quadrature on [0, l] fc via the same change of 
variables (5.2). Further details will be presented elsewhere. 

To see that the orthogonal polynomials tiff, constructed above are eigenfunctions of 
Lj 1 )) 1 , we note that the generalization of (4.9) is 

(5.9) Lfr (x^ 1 ■ ■ ■ = — (\rh\ 2 + (B — 1)|to|^ (x™ 1 ■ ■ ■ x™^ 1 ^ mod 

where B = b\ + ■ ■ ■ + b k +1 and bj = ctj + 1. Normally m k +1 = 0 as the simplex is 
parametrized by any k of the variables x t . Equation (5.9) shows that L^ 1 ™ is upper trian¬ 
gular in any monomial basis ordered by degree (with arbitrary ordering among monomials 
of the same degree), and the eigenvalues can be read off the diagonal. If the Gram-Schmidt 
procedure is applied to these monomials using the dfib inner product, the resulting system 
of orthogonal polynomials will span a nested sequence of invariant subspaces for Lj ™, and 
hence are eigenfunctions, by self-adjointness. The above construction is equivalent to the 
Gram-Schmidt procedure. 


(to > 0), 
(to > 1). 
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Following Griffiths [16], one could combine all the eigenfunctions in each eigenspace 
into reproducing kernels, defined as Gd{x,y) = E|m|=d V’mCaOV’mCy)- In the above ex¬ 
ample, replacing ( x,y,z ) with (xi, X 2 , £3), one would have e.g. Gi{x,y) = 136080 + 
181440[ x\ -x 2 -x 3 -y 1 -y 2 -y 3 + x 1 yi + x 2 y 2 + x 3 y 3 + (zi + x 2 + x 3 )(y 1 + y 2 + y 3 )]. 
This is theoretically appealing due to symmetry and independence of the particular ordering 
chosen to compute the orthogonal polynomials. However, it is more efficient computationally 
to work with the orthogonal polynomial representation directly. It is also numerically unsta¬ 
ble to expand out these polynomials in terms of their coefficients to achieve symmetry in the 
formulas for Gd{x : y) when d is large. 

5.1. Coefficient Estimates. If / is a polynomial, or more generally, a smooth enough 
function on £„, then we can use the recursive approach above to represent / in terms of the 
eigenfunctions of L Klm , {wxP'X/m} ■ Here and in the sequel we assume that all the eigenfunc¬ 
tions are normalized by the condition: 

(5-10) J \il)i,ri,{y)\ 2 dyi, 2 {y) = 1. 

K X 

When I = {j}, then Kx is a vertex; the functions on Kx are constants. We write ipx instead 
of 'f/'x.m, where for all j, x ) = U the constant function. Note that XjV{j} (x) is then a 
function in the null-space of L Klm that is one at e :) and vanishes at all other vertices. We let 
fix. 2 = 5 ej (x) be the atomic measure at the vertex ej of unit mass. 

We begin by subtracting off the values of / at the vertices: 

n+1 

(5.11) f (1) {x) = f(x) - 

j= 1 

The function /l 1 ’ now vanishes at each of the vertices of £ n and can therefore be expanded 
on the 1-skeleton in terms of the eigenfunctions of the form {wx’P'id} ■ where I ranges over 
subsets of {1,... , n + 1} of cardinality 2. The coefficients are computed as integrals over the 
1 -dimensional strata: 

(5.12) c x ^ - [ ^-—^T’ipx,m(x)x(l - x)dx = [ f (1) {x)ipx^{x)dx, 

where x is a coordinate on the stratum Kx- As in (3.12) and (3.14), the rate of decay of these 
coefficients is determined by the smoothness of /l 1 ). 

Subtracting off the contributions from the 1-dimensional strata we get f <2> . which van¬ 
ishes on the 1-skeleton of £„. we can represent this function on the 2-skeleton in terms of 
eigenfunctions of the form {wxipx,rh}, where the cardinality of I is 3. Proceeding in this way 
we get a sequence of functions /, /l 1 ), ..., f^ k \ ..., f^ n \ where /! fe ) vanishes on the 

k — 1-skeleton of 6£„. 

The coefficients coming from f <k) are computed as integrals over the components of the 
fc-dimensional part of the boundary of £ ra : 

15.13) c x ,rn = J ipi:,rh{y)f ( - k) {y)dy. 

Kx 

Here I has cardinality k + 1 and y is a linear coordinate on Kx- Using Shimakura’s for¬ 
mula, (4.1) we can show that 

J L 2 l/Jx,rn(y)f ik) (y)dy = J V’X,m(l/)[Lx™+Kx]/ (fc) (t/)dt/. 

K x Kx 


(5.14) 
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From this we can show that if L 2 ?/ , x,m = Ax.juV’x.m, then, for a function /^ £ C 2l (Kz) 
that vanishes at the boundary of K x , we have: 


(5.15) 


CX.m — 


A' 


V'x.mCy) [Lx ™ +Ki] l f (k) ( y)dy. 


Kx 


From this relation it is clear that the rate of decay of the coefficients {cx,m} is determined by 
the smoothness of f^ on the /c-skeleton, and an L°°-estimate for the eigenfunctions 

To estimate in terms of the original data would take us too far afield, so we conclude 
this discussion by proving sup-norm estimates on the eigenfunctions. The eigenfunctions 
{"fern} are eigenfunctions of the operator Lx. 2 , which is self adjoint with respect to the 
measure d/ix. 2 , and has strictly positive weights. The kernel for the operator e tLx ’ 2 takes 
the formp t (?/, y)dfix^(y), with pt{y, y) = Pt{y , y)- This and the semi-group property easily 
imply that 


(5.16) 


P 2 t(y,y)=[ [ Pt(y,y)} 2 dyx, 2 (y )• 

J K x 


Since the operator has positive weights, we can use the Theorem 5.2 in [8] to conclude that 
there is a constant Ck depending only on the dimension so that 


(5.17) 


P 2 t{y,y) < 


C k 


y-2{B^2i(y)) 


Here B r (y) is the ball in the intrinsic metric (see (2.24) and [8]) of radius y/2t centered at 
y. From the forms of the metric and the measure we can easily show that there are constants 
Co, Ci, so that for small t, on strata of dimension k. we have the bounds 

(5.18) c 0 f fe < MB^y)) < Cifi 
To prove an estimate on ||V’x,m||x°°, we observe that 

(5.19) 


ipx,rh{y) e ' X ’^ = J Pt(y,y)^iMy) d PiAy)- 

Kx 

The Cauchy-Schwarz inequality and the estimates above show that 

C k e~ tXx 


t 2 


(5.20) \^iM\ < 

Setting t = — l/Ax.rn,, gives the estimate 

(5.21) |V’x,m(y)| < C k |Ax,m| 2 • 

Inserting this estimate into (5.15) we see that there is a constant C' k , so that the coeffi¬ 
cients coming from the stratum of dimension k satisfy an estimate of the form: 


lor.ml — C; 


„ ||[L Kim +Kx]'/ (fc) lk" 


(5.22) 
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One-half order of decay is lost with each increase in dimension. We can also give an L 2 - 
estimate, where no such loss explicitly occurs, wherein 


(5.23) 


I Pi,ml < Cl 




A 


i 

X,rh 


These estimates implicitly involve k derivatives of | [L Klm +Kx] l f^ | 2 near the boundary of 
Kx- In both estimates there are further losses that occurs in the estimation of [L Klm +Kx] l f( k ' > 
in terms of the original data /. 

6. The Dirichlet Problem. In many applications of the Kimura equation to problems 
in population genetics one needs to solve a problem of the form 

(6.1) L Klm w = / with u\s= g- 


Here S' is a subset of bE n , generally assumed to be a union of faces. In the constant weight 
case, Shimakura showed that this problem is well-posed if the weights on the faces contained 
in S are all less than 1. Note that, if a weight is greater than or equal to 1, then, with probability 
1, the paths of the associated stochastic process never reach the corresponding face. 

In this section, for simplicity, we continue to consider the case that all weights are zero, 
S = bE n , and that / and g have a certain degree of smoothness. With this assumption, we 
show that the solution to the Dirichlet problem has an asymptotic expansion along boundary 
with the first two terms determined by local calculations, see (6.20) and (6.57). A classical 
example from Population Genetics is the solution to the Dirichlet problem with f = — 1, and 
g = 0, which is given by 

n 

(6.2) u{x) = ^2 ifei "t- 

j =1 0<ii 

where ^(r) = rlogr. For a point x £ E n , the value u(x) is the expected time for a path 
starting at x to reach bE n . This formula is highly suggestive of the form that the general result 
takes; see (6.17) and (6.42) below. 

Boundary data that is not continuous on bE n does arise naturally in problems connected 
with exit probabilities, and cannot be directly treated by the methods in this paper. In these 
types of problems the boundary data is often piecewise constant assuming only the values 
0 and 1, and one has recourse to explicit solutions, see [34]. In this case, one could ap¬ 
proximate the discontinuous boundary data by smooth boundary data, and use the methods 
presented here along with the Feynman-Kac formula in [9] to get upper and lower bounds 
for the solutions of such boundary value problems. For data without some degree of regular¬ 
ity one should not expect the solution to have a simple, explicit asymptotic expansion along 
the boundary. When the solution is not continuous up to the boundary, considerable care 
is required in the interpretation of the partial differential equation along the boundary, and 
boundary condition itself, see [20, 21, 22]. 

To start we consider the simpler problem in a positive orthant S n ,o = R” , with the model 
operator 

n 

(6.3) L 0 = '^2 l x i dl i . 

i=l 

This problem is easier to solve and its solution is nearly adequate to solve the analogous 
problem in a simplex. We start with a simple calculus lemma. 
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LEMMA 6.1. Let ip be a C 2 function of a single variable, then for indices 1 < i\ < ■ ■ ■ < 
ik < n, we have 

(6.4) Loipix^ H-f x ik ) = (x ix H-b x ik )ip"(x^ H-b x ik ). 

The proof is an elementary calculation. 

We now show how to solve the problem 

(6.5) L 0 u = f in S n , 0 with u\ bSn0 = g , 

for / a compactly supported function in C 2 (5 t1j o). We assume that g has a compactly sup¬ 
ported extension, 7j, to S n> o, for which Log also belongs to C 2 (S n ^). It is clear that, generally, 
this problem cannot have a regular solution u £ C 2 (,S' T ,..o). If u belongs to this space, then, for 
any substratum a of dS Ut o, we have 

(6.6) {L 0 u)\ a = L 0 \a u\ a . 

Hence / and g would have to satisfy the very restrictive compatibility conditions 

(6.7) f\a=L 0 \ag\a- 

We look for a solution of the form a = g + v , with v solving 

(6.8) L 0 v = f - L 0 g = f (1) and v \ bSn:0 = 0- 

Our first goal is to give a method for solving (6.8), which gives a precise description of 
the singularities that arise for general smooth, compactly supported data (/, g). We begin 
by giving a precise definition to the meaning of the equations in (6.5): A function u £ 
C 2 (int S n ,o) O C°(S n ,o) solves (6.5) if u\bs n , 0 = and so that L 0 u, which is initially de¬ 
fined only in the interior of the orthant, has a continuous extension to ,5' n .o with Lqu = /, 
throughout. 

In Section 3 we showed how the analogous problem for L Klm is solved in 1-dimension. 
We begin this construction by explaining how to solve (6.5) when n = 1 : For a fixed £<1, 
we let ip £ C£°([0, e)) be a non-negative function, which is 1 in [0, ■§). Suppose that we want 
to solve 

(6.9) xd 2 u = f(x) with u(0) = a 0. 

As noted above, if u £ C 2 ([0,oo)), then, unless /(0) = 0, it cannot simultaneously satisfy 
the differential equation as x —> 0, and the boundary condition. If we write 

(6.10) u(x) = ip{x)[f if))x log x + a] + v(x), 

then 

(6.11) xd 2 u = x-i/> ,, (a;)[/(0)a;loga:-|-a] + 2xip'(x)[f(0){logx + 1)] +ip(x)f( 0) +xd 2 v. 

The first two terms on the right hand side are smooth and compactly supported away from 
x = 0. The function v must solve the equation 

(6.12) xdlv = f(x)-[xip"(x)[f{0)x\ogx + a]+2xip'(x)[f(0)(\ogx + l)]+ip{x)f{0)], 

with u(0) = 0. The right hand side is smooth and vanishes at x = 0, hence the theory 
presented in [6] shows that there is a unique smooth solution to this problem. 
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To carry this out in higher dimensions, we need to introduce some notation. For indices 
1 <«!<•••< ik < n, (or n + 1) we define the vector valued function 

(6-13) -V,,.,, = E x j e j, 

where e 7 are the standard basis vectors for R” (or R” +1 — which is needed will be clear 
from the context). For a function ip defined in R' 1 . the composition satisfies 

( 6 - 14 ) = v\x H =-= x i k = o • 

Our method relies on the following essentially algebraic lemma: 

Lemma 6 . 2 . Let t]{t) = r logr, and let Lq denote the operator defined in ( 6 . 3 ). For 
any k < n, and distinct indices {*i,...,, C {1,..., n + 1}, we have 

( 6 . 15 ) L 0 r ]( xi 1 H- x ik ) = 1 . 


Proof. By the permutation symmetry of the operator Lq it suffices to consider the indices 
{1,..., k}. Since d^r/ir) = 1/t. this follows from Lemma 6.1. □ 

As the next step we write the function v = vq + v\ , where vq will be the “regular” part 
of the solution, vanishing on the boundary, and v\ is the singular part, which also vanishes on 
the boundary and satisfies the equation 

(6.16) £otfit b s„,o=/ ( 1 ) t 6 S„,o ■ 

As we shall see, V\ belongs to C^,p(5 rii o) for any 0 < 7 < 1. In fact we can simply write a 
formula for iq : 

n 

(6.17) = J2 (- i y~ 1 f (1) ( x i 1 ,...,i j M x ii+---+ x i j )- 

j =1 {1 < ii <---< ij < n } 

For example, if n = 2, then 

(6.18) v 1 (xi,x 2 ) = / (1) (aq, 0 ) 77 ( 0 : 2 ) + f (1 \0,x 2 )r]{x 1 ) - / (1) (0, 0 ) 77 ( 2:1 + x 2 ); 

if n = 3, then 

(6.19) 

27(2:1,2:2,2:3) = Z^ 1 (27,2:2,0)77(2:3) + / (1) (xi, 0,2:3)77(2:2) + ,/ (1) (0,2:2,2:3)77(2:1)- 
f (1) (X!, 0,0)77(2:2 +2:3) - / (1) ( 0 , X 2 , 0)77(2:1 +2:3)- 

/ (1) (0,0,2:3)77(2:1 +2:2) + / (1) (0,0,0)77(27 +2:2 +2:3). 


REMARK 6.3. Formula (6.17) should be contrasted to results like those in [4], which 
analyzes boundary value problems for uniformly elliptic operators in domains with comers. 
For this case the classical Dirichlet and Neumann boundary value problems are well posed, 
however, even if the data is infinitely differentiable, these problems typically do not have 
smooth solutions on domains whose boundaries have comer-type singularities. There is no 
analogue of the “regular" solution, which is uniquely determined in the present case. The 
existence of a regular solution is another indication of the remarkable relationship between 
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the degeneracies of the operator L Klm and the singular structure of the boundary of the 
simplex. Singularities arise in the present case from the requirement that the solution satisfy 
a Dirichlet-type boundary condition, which, as we have seen, is generally impossible for a 
function in C 2 (E„). 

A second remarkable feature of this formula is that it provides two locally determined 
terms in the expansion, along 6E„, of the solution u to the original problem in (6.5). We have 

( 6 . 20 ) u(x) = g(x) + vi(x) + 0(dist(a:, 6S„)), 

where it should be noted thatv\{x) = 0(dist(x, 6E n ) logdist(a;, 6E n )). In the non-degenerate 
case, determination of the second term in such an expansion requires the solution of a global 
problem. 

THEOREM 6 . 4 . Let £ C^rp + 7 (E„), then the function zti defined in ( 6 . 17 ) locally 
belongs to C^Z.(S n fi),for any 0 < 7 < 1, as does LqV\. It satisfies the following equations 

( 6 - 21 ) £ 0 ^llbS„,o= f ( 1 ) \bS n ,o and V lt 6 Sn,O =0 • 


Proof. The fact that v\ £ CwfI^.o) follows immediately from the fact that 7 y(r) is 
locally in C^p(S'i i o). We use Lemma 6.2 and the fact that X-^ j. and 27 , + • • ■ + x. tj 
depend on disjoint subsets of the variables to obtain that 
( 6 . 22 ) 


L 0 v 1 = 


3 =1 {l<ii<-*-<i :? <n} 


(-lr 1 Lof^(X, iMxii +---+ Xij ) + fW(X, , ) 


Wt 


It is clear that L 0 v\ £ Cwf(*^«,o) for an y 0 < 7 < 1. The facts that 77 rbs „ i0 = 0, and 
LqVi |"fcs„ 0 = / (1) 0 follow from a direct calculation. We first give the proof for the 

n = 2 case: Observe that 


( 6 . 23 ) v 1 (xi,x 2 ) = / (1) (27,0)77(27) + / (1) ( 0 ,2:2)77(01) - / (1) (0,0)77(27 +2:2). 


We restrict to 27 = 0 to obtain 

(6.24) V! (o:i,0) = / (1) ( 0 , 0 ) 77 ( 2 : 1 ) - / (1) ( 0 , 0 ) 77 ( 2 : 1 ) = 0. 

The case 27 = 0 is identical. Applying Lq we see that 


(6.25) L 0 77(27,27) = xxdzj^ixi, 0 ) 77 ( 2 : 2 )+ 

f { 1 ) (x!, 0) + 2: 2 5 2 2 / (1) (0,27)77(27) + / {1) (0, x 2 ) - / (1) (0, 0). 

Setting 2:2 = 0 now gives 

(6.26) LoVifaO) = f w (xi,0) + / (1) ( 0 , 0 ) - / (1) ( 0 , 0 ) = / (1) (xi, 0 ). 

The general case is not much harder: By symmetry, and continuity it suffices to show 

that 


( 6 . 27 ) 


E E (- i y i f < ' 1) ( X i 1 ,...,i :l )v(xii+---+Xi j )\x n =o=0 and 

3 =1 {l< 2 i<-"<i, <n} 
n 

E E (-» J 1 . + --+ x i i ) + f { 1 Hxi u ...,i i ) 

3 = 1 {l<ii<---<ij<n} 


- 0 


= /(*«)• 
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We prove the first identity in (6.27) by observing that, if ij < n, then X* j. \ Xn — 0 = 
X^i i . We split the sum into two parts and use the fact that r/(x n )\ Xn - o= 0 to obtain: 

n 

(6.28) Y X (~ 1 Y~ 1 f {1 \ x i 1 ,...;i :i )'n(xh + --- + Xi j )\ x „=o= 

3 — 1 {l<ii<-”< 2 j<n} 
n— 1 

j=l {l<ii<"*<ij<n—1} 
n 

X X (- 1 ) i " 1 / (1) ( I i 1 ,..,i J , 1 ,s)i( x i 1 +'-+ x ii- 1 )' 

j —2 {l<ii< -<i 3 _i<n—1} 

Upon changing j — 1 —*■ j in the third line, it becomes clear that the second and third lines 
in (6.28) differ only by a sign, and therefore sum to zero. 

To prove the second identity in (6.27), we observe that the first identity already implies 

that 

n 

(6.29) Y X (-l) i - 1 £o/ (1) (^i 1> ... i i i M*< 1 +"- + * ii )r*„=o=0 J 

3— 1 {l<ii <n} 

so we are only left to prove that 

n 

(6.30) y X (-ir7 (1) (x ?1 , ... A )rx»=o=/(*»)• 

J=1 {l<ii<”-<ij<n} 

Applying the same decomposition as before, we see that 

n 

(6.3D y X (-ir 1 / (1) (7,.,s ] )^=o= 

3= 1 {l<ii<---<ij<n} 

n—1 

/ (1) m + X X L / (1) (^ 1 . hfi)+ 

j =1 {l<ii<---< 2 j<n—1} 
n 

Y X 

j=2 {l<2i<---<z_j_i<n— 1} 

from which the claim is again immediate. □ 

This almost suffices, except for the fact that v\ is not generally compactly supported even 
if the data is. This is because in the definition of Vi we evaluate J 11 > at subsets of the variables 
(cci,..., x„). To repair this we multiply v\ by a specially constructed bump function. To see 
that this works we need another elementary calculus lemma. 

Lemma 6.5. Let ip £ C 2 ([0,oo)) and h £ C 2 (int(«SVi,o))i then 

(6.32) L 0 [hip(xi H-b x n )\ = 

ip(x i H-+ x n )L 0 h + 2tp'(xi H-+ x n )Rh + (xi H-+ x n )-ip"(xi H-+ x n )h. 


Here R = X\d Xl + • • • + x n d Xn . 

Using this lemma we can easily demonstrate the following result: 
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Proposition 6.6. Suppose that /(■*■) e C^p +7 (S'„ j o) is supported in the set {x : xi + 
• • • + x n < N}. If ip € C 2 ([ 0, oo)) equals 1 in [0, N], then 

(6.33) vi(x) = ip{x i H-1- x„)ui(x) 

z.v compactly supported, v\ and LqVi belong to C^p(«S' nj o). The function v\ also satisfies the 
equations: 

(6-34) L 0 vi\bs n ,o= f { 1 ) \bS n , 0 and fi 1 \ bSn0 = 0. 

Proof. The fact that v\ \bs„ 0 = 0 and its regularity properties are immediate. Applying 
Lemma 6.5 we see that 

(6.35) 

L 0 vi = ip(x id- Vx n )L 0 v\+ 2 ip'(xi-\ - \-x n )Rvi + (xi~\ - \-x n )ip" (xH- \-x n )v\. 

The fact that ) = r(logr + 1) implies that Rv i has the same regularity as v\, and 

therefore so does LoVi . Thus to prove the proposition we only need to show that Rv i vanishes 
on the interiors of the hypersurface boundary components. In a neighborhood of the interior 
of Xi = 0 we have the representation v\ = xilogxiai(x) + x±a2(x), where a\ and a2 
belong locally to C 2 . Applying the vector field we see that 

(6.36) Rvi = a:i(loga;i + l)ai(x) + x± \ogx\Ra\{x) + xi(a,2(x) + RCL2), 

which obviously vanishes as X\ —> 0. The other boundary faces follow by an essentially 
identical argument. □ 

We now can complete the solution of (6.8) and (6.5) as well. Write v = vq + vi, where 

(6.37) L 0 vo = f (1) - L 0 vi and f 0 \ bSn O =0. 

Since /« - LqVi g (S„, 0 ) is compactly supported and vanishes on the boundary, it 
follows from the results in [7] that (6.37) has a unique solution vq G C^,p +7 (5 raj o)- While 
we call vo the regular part of the solution, it will not in general be smooth. It will have 
complicated singularities of the form [x^ H— • + Xi k \ l [log(xi 1 +■ —h Xi k )] m , for 2 < l, and 
with the maximum value of to bounded by a function of l. 

With small adaptations this method can also be used to treat the case of the n-simplex, 
and the neutral Kimura diffusion operator, L Klm . As noted above, the n-simplex is most 
symmetrically viewed in the affine plane xi + ■ ■ ■ + x n +i = 1. This representation makes 
clear that every vertex is identical to every other vertex. For the construction of the solution 
to the Dirichlet problem: 

(6.38) L Klm u = / in £ n with u\bs n = g 

it turns out to be simplest to work in the non-symmetric representation, with one vertex iden¬ 
tified with 0. Recall that 

(6.39) £„ = {(xi, ..., x n ) : 0 < x», i = 1,..., n ; x\ 4-+ x n < 1}. 

It is clear that there is a choice of linear projection into R" so that any given vertex of £„ is 
so identified with 0. For the moment we work in £„ with coordinates (xi,..., x n ). 

We begin by assuming that the data /, g is supported in a set of the form = {x : 

0<Xi+ - —hx n <l — 2e}, for some 0 < e. As before we let g denote an optimally smooth 
extension of g as a function with support in £ n>e . We write u = g + v, where v satisfies 

L Kim u = /-L Kim 5 = /W with xr h£n =0. 


(6.40) 
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As before we write v = vo + v ±, where 

(6.41) vi = i/j(x i H-h x n )vi(x), 

with %j) G C^°([0,1 — |)), satisfying tp( T ) = 1) for T £ [0,1 — e]. The use of the bump 
function is much more critical here, as we do not have good control on the function v\ on the 
face where x\ + • • • + x n = 1. 

The function v\ is defined by the sum: 

n 

(6.42) ^iO)=E E (-1) J_1 / (1) (^ 1 , + ''' + 

3 =1 {l<ii<---<ij<n} 

What is special about the choice of coordinates is the possibility of having the two functions, 
and V&ii + • • • + Xij ), depend on disjoint sets of coordinates. Elementary 
calculations show that 

(6.43) L Klm t]{x ix 4-f x i:) ) = 1 - Oj, 4-f x i} ), 

and 


(6.44) 

L Kim^(i)( Xji + ... 4 -ajfJ] = rj(x il +■■■ + x ij )L Kim f w (X lu j.)- 

2( Xil + ■ ■ ■ + Xi^rffa + ■ ■ ■ + .r' s-V,^ ^) + 

4- (1 - {Xi ± 4-f X; h ^), 


where the vector fields R~- $ are defined by 


(6-45) -4 = E 

&£{*<> — ,41 

Let = bYi n \ {a: : x\ 4- ■ • • 4- x n = 1}. Arguing as before we show that v\ satisfies 
the equation along bYf n . 

Theorem 6.7. Let /W G C^p +7 (E n ). The function v\ defined in (6.42) belongs to 
C^rp(E n ), as does L Klm v\. It satisfies the following equations 

(6.46) LKim ^r 6 £ ; =/ ( 1 ) r 6 g ; and fi[ 6 g^=0. 


Proof. The regularity statements for v\ and L Klm v\ follow easily from (6.42) and (6.44). 
The fact that v\ f b g, = 0 follows from Theorem 6.4. The proof that L Klm vi [ b g, = /W [ b g, 
is similar to the proof of Theorem 6.4. As before continuity shows that we only need to prove 
this statement for the interiors of the hypersurface faces, and symmetry shows that it suffices 
to consider {x n = 0}. The terms coming from the first line of (6.44) can be treated exactly 
as before. For the terms from the third line of (6.44), the sole difference is the coefficient 
(1 — x n ), which equals 1 where x n = 0. The first order cross terms, coming from the second 
line of (6.44), require some additional consideration. 

We once again use the observation that if ij < n, then ;• [ Tri=0 = X-^ 5 . as 

well as the facts that 


%.s/'At,.;,)(—»= 




,n/ 
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and 


\x n — 0— {x n X n T X n ) \x n — 0— 0; 

to split the sum into two parts, obtaining: 


(6.47) 

n 

Y Y (- 1 ) J_l2 ( ;r ii+- • ■+x ij )r}'{ Xil +- ■ ■+x ij )Ri 1 ^ jJ (1 \X iu j j )\ Xn= o 

j= 1 {l<u<---<z :? <n} 
n—1 

H H (-l) J_1 2(a; il +- • •+a: i ,)r? , (a; il +- • ■+x ij )R- i ^j.~f il \X ii }i . fi )+ 

j= 1 {l<ii<"-<ij<n—1} 
n 

Y Y (- 1) J ' _1 2 (Si!+- • •+^-)t/(x il +- • 

J—2 {l<ii<--.<ij_i<n—1} 


As before, after changing variables in the third line with j K > j — 1, we see that the second 
and third lines differ only by a sign, demonstrating that these terms sum to zero along the face 
given by {x n = 0}. This completes the proof of the theorem. □ 

To complete the discussion of this case we need to check that Vi also satisfies the equa¬ 
tions in (6.46). An elementary calculation shows that 


(6.48) L Klm [uiV’(:ci H-ha; n )] = 

4>(xi 4-+ x n ) L Klm v\ + H-1- x n )Rv i + vi L Klm i/j, 


where 

n 

(6.49) R = 2(1 — (x\ 4 - \-x n ))'^ j x i d Xi . 

7=i 

As before, R is tangent to 6E„, and Rv i £ C^p(E„). It is easy to see that 

(6.50) Rvi\ b ^, = 0. 

Summarizing, we have shown 

PROPOSITION 6.8. Let /(■*■) £ C^rp +7 (S ra ). The function v\ belongs to C^rp(E ra ), as 
does L Klm ui. Iff)(x\ + ■ ■ ■ + Xn ) = 1 on supp C E n ^,for an e > 0, then vi satisfies 
the following equations 

(6.5!) LKlm ^r b£n =/ (1) r fe£n and 


We can now complete the solution of (6.40). We write v = vo + v\, the function vo must 
satisfy 

(6.52) L Kim uo = / (1) — L Kim ui with t)o|'bs„= 0. 

The existence of a solution vo £ C^,p +7 (E„), which can be taken to vanish on the boundary, 
follows from the results in [7], Numerically, the method of Section 5 (with k = n) can be 
used to solve (6.52) for vq. 
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Finally if we have general data (f,g) so that /, and L Klm gi £ C^p +7 (E„), then we 
choose a partition of unity {pi,..., ip n+ 1 } so that the function <pj equals 1 in a neighborhood 
of the vertex ej C £ n C R" +1 , and vanishes in a neighborhood of the opposite face, where 
Xj = 0. The data {pj(f, g) : j = 1,... ,n + 1} satisfies the hypotheses used above, and 
therefore we can construct solutions {uj} to the equations 

(6.53) L Kim Uj = ipjf on £„, with uj\b'£ n = Pjg, for j = 1,..., n + 1. 

To construct the solution in a neighborhood of ej £ R n+ 1 we project the simplex into the 
hyperplane {Xj = 0}, and use the projective representation where e,j corresponds to 0. 

Now setting 

(6.54) u = ui H-+ zz n +i, 

we obtain a solution to the original boundary value problem 

(6.55) L Klm zt = / on £„ with z4b£ n = g. 

This proves the following general result: 

Theorem 6.9. Let f £ C^rp +7 (E„), and g have an extension g to £ n so that L Klm g £ 
C^rp +7 (E„). The Dirichlet problem (6.55) has unique solution u, which takes the form u = 
g + u + z/ 1 ), where £ C^rp +7 (£ n ), vanishes on the boundary. The singular part, u ^ 
is of the form 

n+1 n 

(6.56) zt^(x) = ££ E . +---+X.,) 

i= 1 j =1 {l<zi<n+l: 

Here F is a function constructed from the pullbacks of {'Pj f^ 1 } to the affine model £„. 

REMARK 6.10. The solution u to the boundary value problem in (6.55) also has an 
explicit 2-term expansion at the boundary similar to that given in (6.20): 

(6.57) u(x) = g(x) + zz^^(x) + 0(dist(x, bH n )). 


Proof. Everything has been proved except the uniqueness statement. This follows from 
the maximum principle and the facts that u is continuous in the closed simplex, and L Klm is 
a strongly elliptic operator in the interior of the simplex. □ 

We observe that u has, in some sense, very complicated singularities, in that it includes 
a smooth function times (x^ + ■ ■ ■ + Xj ,) log(xi i + • • • + Xi .), for each set of indices 

1 < ii < ■ ■ ■ < ij < n + 1 , for j £ {1, ..., n}. 

To resolve these singularities to be classically conormal, one would, in principle, need to 
successively blow up all the strata of the boundary, starting with the codimension n parts 
and proceeding upwards to the codimension 2 part. Given the very explicit form that this 
singularity takes, such an approach would only obscure its simple and rather benign structure. 
This general approach is pursued in the papers [20, 21, 22]. These authors do not require the 
data (/, g) to be continuous, but allow data with complicated singularities along the boundary. 

We remark that the approach proposed in Section 4 for solving the regular problem 
L Klm u = f (without boundary conditions) can be modified slightly to provide an algorithm 
for extending g (defined on the boundary) to g (in the interior) in such a way that L Klm g is 
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very easily computed. We write g = go + • • • + g n _ 1? where go = Y^=i 9( e j) x j a S rees with 
g on and 

k -1 

i=o 

The right-hand side is zero on E„ - 1 , so this equation decouples into independent homoge¬ 
neous “Dirichlet” extension problems from the connected components of E „ \ E^ _1 to E n . 

Using the eigenfunction expansion to represent the right-hand side on each face of Ejj 
gives the desired representations 

(6.58) gk = and L * g ^ 

I I 

with both functions canonically defined throughout E n . If g is a polynomial of degree less 
than or equal to d on each face of 6E n , then g £ Vd as well. Thus, we can replace / by 
/(i) = / - L Kim g in (6.40) at the outset, thereby avoiding non-homogeneous boundary 
conditions when working with the partition of unity. 

Our approach to solving the Dirichlet problem works equally well if we add a vector field 
V that is everywhere tangent to 6E„. In a projective chart, such a vector field takes the form 

n 

(6.59) V = {x)xjd Xj , 

3 =1 


with the additional requirement that 

n 

(6.60) ^ x :i ( x ) [ x i h—— 0 - 

i=i 


A simple calculation shows that 


n 

(6.61) Vriix^ H-b x id ) = ^ ai{x)xi A {ilj ... )ij .}(()[log(x il H-fa:^) + 1 ], 

z=i 


where 


(6.62) 




lif l G {ii,... ,ij} 

0 if l £ 


The function on the right hand side of (6.61) is easily seen to be continuous on the closed 
n-simplex. In fact these functions belong to C^p(E n ), for any 0 < 7 < 1. Applying V to vi 
we see that 


(6.63) Vvi = viVip(xi H-b x n ) + H--b x n )Vv\, 


from which it is clear that Vv± is continuous on S„. To show that Vv\ t&E„ = 0, it there¬ 
fore suffices to prove it in the interiors of the hypersurface boundary faces, but this is clear 
from (6.61) and (6.63). With these observations it follows that (L Klm -bU)t>i £ C^p(E n ) 
and 


(6.64) 


(L Kim +V)V! = L Kim V! \ bSn = / (1) 
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Proceeding as above we easily demonstrate that, if V is tangent to the boundary b£ n , 
then the Dirichlet problem: 

(6.65) (L Kim +V)u = / in E„ with u\bx n = g 

has a unique solution of the form u = Uq + u±, where ui is given by the formula in (6.56), 

and u 0 e C^p +7 (£„). 
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